A general purpose Fortran 90 electronic 
structure program for conjugated systems using 
Pariser-Parr-Pople model 



Priya Sony\ Alok Shiikla^ 

Department of Physics, Indian Institute of Technology, Bombay, Powai, Mumbai 

400076, INDIA 



Abstract 

Pariser-Parr-Pople (P-P-P) model Hamiltonian has been used extensively over the 
years to perform calculations of electronic structure and optical properties of vr- 
conjugated systems successfully. In spite of tremendous successes of ab initio theory 
of electronic structure of large systems, the P-P-P model continues to be a popular 
one because of a recent resurgence in interest in the physics of vr-conjugated polymers, 
fullerenes and other carbon based materials. In this paper, we describe a Fortran 90 
computer program developed by us, which uses P-P-P model Hamiltonian to not only 
solve Hartree-Fock (HF) equation for closed- and open-shell systems, but also for 
performing correlation calculations at the level of single configuration interactions 
(SCI) for molecular systems. Moreover, the code is capable of computing linear 
optical absorption spectrum at various levels, such as, tight binding (TB) Hiickel 
model, HF, SCI, and also of calculating the band structure using the Hiickel model. 
The code also allows the user to solve the HF equation in the presence of finite 
external electric field, thus, permitting calculations of quantities such as static 
polarizabilities and electro-absorption spectra. Additionally, it can perform transformation 
of P-P-P model Hamiltonian from the atomic orbital (AO) representation (also called 
site representation) to the molecular orbital (MO) one, so that the transformed 
matrix elements can be used for high level post-HF calculations, such as, full CI 
(FCI), quadruple CI (QCI), and multi-reference singles-doubles CI (MRSDCI). We 
demonstrate the capabilities of our code by performing calculations of various properties 
on conjugated systems such as frans-polyacetylene (t-PA), poly-para-phenylene (PPP), 
poly-para-phenylene-vinylene (PPV), 0%0-acenes, and graphene nanodisks. 
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Program Summary 

Title of program: ppp.x 
Catalogue Identifier: 
Program summary URL: 

Program obtainable from: CPC Program Library, Queen's University of Belfast, 
N. Ireland 

Distribution format: tar.gz 
Computers : PC's/Linux 

Linux Distribution: Code was developed and tested on various recent versions 
of Fedora including Fedora 11 (kernel version 2.6.29.4-167) 
Programming language used: Fortran 90 

Compilers used: Program has been tested with Intel Fortran Compiler (non- 
commercial version 11.1) and gfortran compiler (gcc version 4.4.0) with opti- 
mization option -O. 

Libraries needed: This program needs to link with LAPACK/BLAS libraries 
compiled with the same compiler as the program. For the Intel Fortran Com- 
piler we used the ACML library version 4.3.0, while for gfortran compiler we 
used the libraries supplied with the Fedora distribution. 
Number of bytes in distributed program, including test data, etc.: size of the 
tar file bytes 

Number of lines in distributed program, including test data, etc. : lines in the 
tar file 

Card punching code: ASCII 

Nature of physical problem: The problem of interest at hand is the electronic 
structure of vr-conjugated systems. For such systems, the effective 7r-electron 
P-P-P semi-empirical model Hamiltonian proposed by Pariser, Parr, and Pople 
offers an attractive alternative as compared to the ab initio approaches. The 
present program can solve the HF equations for both open- and closed-shell 
systems within the P-P-P model. Moreover, it can also include electron cor- 
relation effects at the singles CI level. Along with the wave functions and 
energies, various properties such as linear absorption spectra can also be com- 
puted. 

Method of Solution: The single-particle HF orbitals of a 7r-conjugated system 
are expressed as linear combinations of the p^-orbitals of individual atoms 
(assuming that the system is in the xy-p\ane). Then using the hopping and 
Coulomb parameters prescribed for the P-P-P method, the HF integro-differential 
equations are transformed into a matrix eigenvalue problem. Thereby, its solu- 
tions are obtained in a self-consistent manner, using the iterative diagonalizing 
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technique. The HF orbitals thus obtained can be used to perform a variety of 
calculations such as the SCI, linear optical absorption spectrum, polarizabilty, 
electro-absorption spectrum, etc. 

Running Time: The examples provided each only take a few seconds to run. 
For a large molecule or a polymer, however, the run time may be a few min- 
utes. 

Unusual features of the program: None 



1 Introduction 

The studies of electronic structure and optical properties of 7r-conjugated sys- 
tems have fascinated physicists and chemists alike for a long time[l,2], be- 
cause of their possible applications in opto-electronic devices such as OLEDs, 
OFETs, organic lasers, solar cells, e^c[3]. Ever since graphene and its nano- 
structures such as graphene nano-ribbons, nano-disks, etc. have been synthesized[4], 
interest in the study of 7r-conjugated systems has multiplied many folds[5]. 

In the 1950's, Pariser, Parr, and Pople proposed a simple model which incor- 
porates the essential physics of correlated conjugated systems in a very elegant 
manner [6]. In this model only vr-electrons are considered explicitly, while the 
effect of (j-electrons is included in an implicit manner in terms of various pa- 
rameters. Moreover, long range electron-electron interactions are taken into 
account by means of suitable Coulomb parameters. Thus, one reduces the 
number of degrees of freedom involved in the underlying Schrodinger equa- 
tion considerably. Because of the lack of large scale computational facilities 
during the 1950's such an approach was unavoidable even for relatively small 
molecules such as benzene. However, the remarkable fact is that in spite of 
so many approximations involved, P-P-P model based calculations were ex- 
tremely successful in describing the electronic structure, in general, and the 
optical properties of such systems, in particular [2]. 

Now, the question arises that in the present times, when computational power 
has grown by orders of magnitude as compared to those in 1950's, why should 
one be still interested in using the P-P-P model, particularly, when so many 
sophisticated ab initio electronic structure programs are available. In our opin- 
ion there are two main reasons behind the continued use of the P-P-P model 
at present: (a) One gets to understand the electronic and optical properties 
of such systems within a rather simplified picture, in terms of the dynamics 
of the TT-electrons, (b) Ever since the revival of interest in the 7r-conjugated 
polymers because of their possible applications in opto-electronics, the size 
of the systems being investigated these days is becoming quite large consist- 
ing of hundreds of atoms, pendant side groups, and terminal functionalities. 
Therefore, simultaneous high quality ab initio description of both the ground 
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and the excited states of such systems (needed to describe their hnear and 
nonlinear optical properties) is virtually impossible. While, the P-P-P model, 
with its reduced degrees of freedom, allows one to compute low-lying excited 
states of rather large systems using the present computational capacities. 

Indeed a large number of calculations on various conjugated systems have 
been performed in last several years using the P-P-P model Hamiltonian 
with a great deal of success[l,7,8,9,10,ll,12,13]. In our own group and with 
our collaborators, we have performed several calculations on systems as di- 
verse as phenyl disubstituted polyacetylenes (PDPAs)[14,15,16,17,18,19,20], 
PPVs[21,22,23], and o%o-acenes[24,25,26]. In all the afore-mentioned works 
the HF calculations and the transformation of the Hamiltonian from the AO 
to the MO representation were performed using a Fortran 77 code, developed 
earlier in our group[27]. However, the code in question lacked the capability 
of performing unrestricted Hartree-Fock (UHF) calculations, and thus, was 
limited only to closed-shell systems. Moreover, Fortran 77 language did not 
allow dynamic memory allocation, resulting in a dependence on compile-time 
dimensional parameters. Therefore, we have rewritten the code in a modern 
programming language, namely, Fortran 90, which allows it to utilize dynamic 
memory allocation, thereby freeing it from various array limits, and resultant 
artificial restrictions on the size of the molecules. Thus our program can be 
used on a given computer until all its available memory is exhausted. The 
present computer program can perform restricted Hartree-Fock (RHF) calcu- 
lations on closed-shell systems, as well as unrestricted-Hartree-Fock (UHF) 
calculations on open-shell systems. Moreover, like the earlier version of the 
code[27], it can also carry out correlation calculations at the level of singles 
configuration interactions (SCI) for molecular systems. Additionally, the code 
is capable of computing optical spectrum at various levels, such as, TB, HF, 
and SCI, and also of calculating the band structure using Hiickel model. The 
program also allows user to solve HF equation in the presence of an exter- 
nal electric field, thus, permitting calculation of dielectric polarizabilities and 
electro-absorption spectrum of conjugated molecules. Additionally, it can per- 
form transformation of P-P-P model Hamiltonian from the AO representation 
(site representation) to the MO one, so that the results can be used for high 
level post HF correlated calculations, such as, full CI (FCI), quadruple CI 
(QCI), and multi-reference singles-doubles CI (MRSDCI). The other features 
of the code include the charge density analysis and the molecular orbital anal- 
ysis. Apart from describing the computer program, we also present several of 
its applications which include various examples for t-PA, PPP, PPV, acenes, 
and graphene nanodisks. 

The remainder of the paper is organized as follows. In section 2 we briefly 
review the theory associated with the P-P-P model Hamiltonian. Xext, in 
section 3 we discuss the general structure of our computer program, and also 
describe its constituent subroutines. In section 4 we briefly describe how to 
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install the program on a given computer system, and to prepare the input 
files. Results of various example calculations using our program are presented 
and discussed in section 5. Finally, in section 6, we present our conclusions, 
as well as discuss possible future directions. 



2 Theory 



In this section we briefly review the theory associated with the P-P-P model 
Hamiltonian and provide the necessary HF equations within the basis set 
approach. 



2.1 Hamiltonian 



The P-P-P Hamiltonian can be written as 



Hppp = J2 (^i4aCia + ^iji4aCja + cJ^Q^,) + U ^i^^il 

i,<T {ij),(7 i 

+ T.Vvini-l){nj-l), (1) 

i<j 

where represents the site energy associated with the ith atom, (ij) implies 
nearest neighbors, cj^ creates an electron of spin a on the orbital of carbon 
atom i (assuming that the system is lying in the xy-plane), rii^ = cj^Cicr is 
the number of electrons with spin a, and rii = TifjUia- is the total number 
of electrons on atom i. The parameters U and Vij are the on-site and long- 
range Coulomb interactions, respectively, while Uj is the nearest-neighbor one- 
electron hopping matrix element. 

The particular forms of the parameters U and Vij determine the specific vari- 
ants of the P-P-P model Hamiltonian[7]. For U 0, and Vij = 0, the model 
becomes the Hubbard Hamiltonian, while for U 0, and Vi 7^ 0, but all 
Vij — 0, we obtain the extended Hubbard model. Choosing for the long-range 
Vij the form 

_ U 
^ 7 \ . 2-1 1/2 ' 



with U and ro fitted to data gives the Ohno variant[28] of the P-P-P model, 
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whereas taking 



u 



(3) 



gives the Mataga-Nishimoto variant [29]. In exponential variant [30], Vij takes 
the following form 

V^. = ^exp(-^) , (4) 



In above Eqs.(2, 3, 4), Rij — jr, — Vj\ is the distance between sites i and j in 
A, while To is some parameter in the same units. 

In this work, we report calculations based upon the Ohno parametrization of 
the P-P-P model mentioned above (c/. Eq.2). Moreover, to account for the 
inter-chain screening effects we use a slightly modified form introduced by 
Chandross and Mazumdar[31], 

y,,- = C//«;,,(l + 0.61174)'/' , (5) 



where Kij depicts the dielectric constant of the system which can simulate the 
effects of screening and Rij is defined above. In various calculations performed 
on phenylene-based conjugated polymers including PDPAs[14,15,16,17,18,19,22, 
it was noticed that "screened parameters" with U = 8.0 eV and Ku — 1.0, and 
Kij = 2.0, otherwise, proposed by Chandross and Mazumdar[31], lead to much 
better agreement with the experiments, as compared to the "standard param- 
eters" with U = 11.13 eV and Hij = 1.0, proposed originally by Ohno. In the 
present code we provide the user with the freedom to choose these "standard", 
"screened", or any other user defined parameters for the Coulomb interactions. 

In order to calculate static polarizabilities, one can solve the HF equations in 
the presence of an external static electric field. Thus, we have modified Eq.(l) 
under the electric dipole approximation by introducing the corresponding term 
containing the uniform electric field E. The overall Hamiltonian of the system 
is then given by 

H^J^f = Hppp - fi.E = Hppp + |e|E.r , (6) 



where Hppp is the unperturbed Hamiltonian that describes the system in the 
absence of the electric field, e represents the electronic charge, n = — er, is 
the dipole operator, and r is the positive operator. 
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2.2 UHF equations 



Now we describe our formulation for the UHF method. The corresponding 
RHF equations can be easily deduced from them. As per the assumptions of 
the UHF method, we assume that the ;U-th up- and down-spin orbitals are 
distinct, and are represented (say) as and ip'^^\ respectively. We assume 
that these orbitals can be written as a linear combination of a finite-basis set 

<"^ = E0^' (7) 



where 0i's represent the p^-orbitals in question, and the determination of the 
unknown coefficients C^^ amounts to the solution of the UHF equations. It 
is further assumed that the 0j orbitals form an orthonormal basis set. In the 
equation above, we have only stated the expressions for the up-spin orbitals, 
the case of the down-spin orbitals can be easily deduced. We further assume 
that the total number of up-/down-spin electrons is N^/Np, such that Na + 
Np — Ne. Using the conjecture of Eq. (7) in conjunction with the Hamiltonian 
above, one obtains the so-called Pople-Nesbet equations 

E(^" - = 0, (8) 

j 



where is the UHF eigenvalue of the /x-th up-spin orbital, is the Fock 
matrix for the up-spin electrons defined by the equations 

F^ = U,-P^V,,, (t^j) (9) 

FS-e,-Y.yij+T.PjjVij-PSVu, (i^j) (10) 



where e^, tij and Vij are defined above (c/. Eq. 1), and and Pjj, are the 
up-spin and total density matrix elements, respectively, defined as 

No, 

p^ = T.cff*c^^ (11) 



and 



p. . = P." + PK 



:i2) 



Once the Fock matrix is constructed, one solves the eigenvalue problem for the 
up-spin Fock matrix (c/. Eq. 8) as well as the down-spin Fock matrix, using 
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the iterative diagonalization technique, to achieve self-consistency. From the 
equations given above, it is easy to deduce the expressions for Ffj^ as well as 
the Fock matrix elements for the RHF case. 



3 Description of the Program 



Our computer code consists of the main program, and various subroutines, all 
of which have been written in Fortran 90 language. Additionally, the program 
must link to the LAPACK/BLAS library, whose diagonalization routines are 
used by our program. 

Before we describe various subroutines in detail, it is fruitful to discuss the 
scaling of the memory requirements, and the cpu time (calculation time), with 
the system size. If N is the number of sites (tt orbitals) in the system, then 
the memory requirements will roughly scale as N"^, needed mainly for the 
storage of the Fock matrix, and its eigenvectors. The cputime, on the other 
hand, will be dominated by the self-consistency iterations (the so-called "self- 
consistent field" (SCF) process), which consists of construction of the Fock 
matrix, and its diagonalization. The cputime needed for the construction of 
the Fock matrix scales as N"^, while the diagonalization time exhibits 
scaling. Thus, we conclude that the dominant scaling of the cputime is N^, 
with respect to the system size. For the open-shell systems, for which the UHF 
procedure is needed, the memory and cpu time requirements will roughly be 
twice as compared to a closed shell system of similar size because the UHF 
procedure involves the construction and diagonalization of two Fock matrices, 
corresponding to the "up", and the "down" spins. 

Several molecules in nature are highly symmetric, and, therefore, in general, 
it is fruitful to implement the point-group symmetries in electronic-structure 
codes like the present one. However, our code does not make use of point- 
group symmetries mainly because the computation time associated even for 
large systems is really insignificant. Therefore, the programming effort re- 
quired to implement various irreducible representations will simply outweigh 
the rewards in terms of reduced cpu time. Moreover, the irreducible repre- 
sentations associated with various orbitals and many-particle states can be 
easily determined by examining the transition dipole moments among them, 
in conjunction with the dipole-selection rules of various point groups. Next, 
we briefly describe the main program, and each of the subroutines of our code. 
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3.1 Main program PPP 

This is the main program of our package which reads input data such as 
which Hamiltonian to use, its parametrization, charge on the system, total 
number of atoms in the unit cell, and their coordinates. The program also 
reads in the options to perform various types of calculations such as RHF, 
UHF, SCI, optics, etc. The option is also provided to perform HF calculations 
in the presence of an external electric field, which is specified by its Cartesian 
components in the units of V/A. The program can delete the atoms and cal- 
culates the entities like total number of sites, electrons and occupied orbitals. 
It dynamically allocates various arrays, and then calls other subroutines to 
accomplish the remainder of the calculations. Because of the dynamical array 
allocation, the user need not worry about various array sizes, as the program 
will automatically terminate when it exhausts all the available memory on the 
computer. 

3.2 Subroutine R_ ATOM 

This subroutine reads the coordinates of the atoms of a unit cell with respect 
to its user specified origin. It can also generate coordinates of some important 
structural units such as a bond, phenyl group (benzene), and even fullerene 
(Ceo) to facilitate an easy realization of the molecular system under consider- 
ation. 

3.3 Subroutine BENPERP 

This routine generates the coordinates of six carbon atoms forming the ben- 
zene backbone. The orientation is perpendicular to the conventional orienta- 
tion. Note that the origin of the benzene skeleton is considered to be the center 
of the hexagonal ring. Rotations, if desired, are performed keeping the origin 
fixed. Positive angles are considered to be anti-clockwise, while the negative 
angles are treated clockwise. Finally, the center is translated to the location 
specified by the user. The plane in which benzene locates is taken as a card 
(XY or YZ or ZX) in the input file. 

3.4 Subroutine BEN ZEN 

Similar to subroutine BENPERP, this routine also generates the coordinates of 
six carbon atoms forming the benzene backbone. This routine differ from the 
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previous one in a way that it creates benzene ring in the orientation parallel 
to the conventional orientation. One can also rotate and translate the ring, 
and can also specify the plane in which benzene ring lies. 



3.5 Subroutine BOND 



This routine generates a bond (2 carbon atoms). Center of the bond is taken 
to be the origin and initially the bond is assumed to be along the x-axis. 
Rotations, if desired, are performed keeping the origin fixed. Positive angles 
are considered for anti-clockwise rotation, while the negative angles are treated 
for clockwise rotation. Finally, the center is translated to the location specified 
by the user. 



3.6 Subroutine STLINE 



This subroutine also generates a bond, but with the first atom of the bond at 
the origin. Here, initially the bond is assumed to be along the x-axis. Opera- 
tions like rotation and translation can also be performed. 



3.7 Subroutine C60 GEN 



This routine, which is originally written by Dharamvir et a/. [32], generates the 
position coordinates of all 60 carbon atoms of the fullerene structure. Inputs 
are two bond lengths, i.e., single and double bond lengths. This does not have 
the "standard" orientation, but one with pentagons at top and bottom. 



3.8 Subroutine R SITE 



If the molecule under consideration consists of translationally invariant units, 
this routine generates the coordinates of all the 7r-electron sites in the system 
from the coordinates of origins of the reference unit cell and the translational 
vector of the lattice. At present this subroutine is restricted to periodicity only 
in one dimension, thus, making it useful for polymers and their oligomers. 
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3. 9 Subroutine DEL A TM 



This subroutine deletes the users specified atoms from the list containing all 
the atomic coordinates. 



3.10 Subroutine SYMSITE 

This routine symmetries the coordinates of the atoms to place the origin of 
the molecule at the center of mass of the molecule. 



3.11 Subroutine PRINTR 

This routine truncates the position coordinates to seven places to right of the 
decimal and then prints out the coordinates of all the sites into the output 
file. In addition, it also creates a file named as 'atomic_coord.xsf', which can 
be used as input file in XCrySDen[33] program, meant for visualizing the 
molecular structure of the system under consideration. 



3.12 Subroutine MATEL 

This is the master routine meant for generating the one-electron (ty) and 
two-electron (V^j) matrix elements, with or without the electric field. This is 
done based upon the data specified by the user in the input file such as the 
Hamiltonian under consideration. Coulomb parameters to be used (if any), 
hopping matrix elements connecting various sites, etc. 



3.13 Subroutine Tij READB 

This routine is meant for reading the hopping related input for the band 
structure calculations using the TB model (Hiickel model). It reads the unique 
intracell, as well as intercell, hopping matrix elements, and their connectivities 
from the input file. Note that an one dimensional tight-binding system with 
only nearest-neighbor coupling is assumed. Thus, a given cell can only connect 
to at the most to two neighbors (one to the left and the other to the right). 
For anything else the code will have to be modified. Also one has to specify 
the unit cell to which it is connected. 
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3.14 Subroutine READ _EI 



If there any site energies are, this routine reads them. In case of more than one 
unit cell, the routine also generates rest of the site energies using translational 
invariance. 



3.15 Subroutme Tij_READ 



This subroutine reads the unique hopping matrix elements and their connec- 
tivities from the input file. In case the system consists of more than one unit 
cell, the routine also generates rest of the hoppings from translational invari- 
ance. 



3.16 Subroutine Tij_ GEN 



This routine is used in order to generate hopping automatically from the 
information of bond distance and the corresponding hopping matrix elements. 
The subroutine considers all pairs of sites in the molecule and if the distance 
between them is equal to the user provided distance, it assigns to them the 
user provided hopping matrix elements. This automates the generation of 
hopping matrix elements to a great extent, thus, simplifying the input data. 
This routine is invoked by the card 'HOPGEN' in the input file. 



3.17 Subroutine Vij_CAL 



This subroutine computes the long-range Coulomb matrix elements Vij {cf. 
Eqs. 2, 3, & 4) for the P-P-P model with three possible parametrization namely 
Ohno, Mataga-Nishimoto, or exponential, depending upon the input provided. 
In case the choice of Hamiltonian is the Hubbard or the extended Hubbard 
model, the relevant Coulomb parameters U and/or V, are provided in the 
input itself. In case of RHF calculations, the routine can also calculate the 
contribution corresponding to the correlation functions, if such a calculation 
has been requested in the input file. For any choice of a correlated Hamiltonian, 
corresponding interaction matrix elements are stored in an array. 
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3.18 Subroutine BANDS 



This is the subroutine which performs the band structure calculations using 
the Hiickel model. In order to get the band structure, first Fourier transforma- 
tion of the hopping matrix elements is performed, and later the matrix is diag- 
onalized using the subroutine ZHPEV from the LAPACK/BLAS library to ob- 
tain the eigenvalues and eigenvectors. The eigenvalues of corresponding bands 
are written in an ASCII file 'bands.dat' which can be used for plotting the 
bands using standard programs such as xmgrace[34] or gnuplot[35]. The eigen- 
vectors are written in the binary format, in the file named 'bloch_orbitals.dat'. 



3.19 SubrouUne FOUTRA TB 



This subroutine computes the Fourier transform of one dimensional nearest- 
neighbor-interacting tight-binding Hamiltonian. It is called by subroutine BANDS, 
discussed above. 



3. 20 Subroutine SCF_ RHF 



This subroutine solves the RHF equations for the system under consideration 
in a self-consistent manner, using the iterative diagonalization procedure and 
returns the canonical SCF orbitals, their eigenvalues, and the total energies. 
The arrays which are needed during the calculations are allocated before the 
calculations begins, and are deallocated upon completion. Before the first it- 
eration, Hiickel model Hamiltonian is diagonalized to obtain a set of starting 
orbitals. Subsequently, the Fock matrix corresponding to those orbitals is con- 
structed and diagonalized. The process is repeated until the self-consistency is 
achieved. During the self-consistency iterations, subroutine DSPEVX from the 
LAPACK/BLAS library is used to obtain the occupied eigenvalues and eigen- 
vectors. Obtaining only the occupied eigenpairs, as against the entire spec- 
trum, leads to considerable savings of CPU time for large systems. However, 
if the entire spectrum of eigenvalues and eigenvectors is needed, say, to perform 
optical absorption calculations or preparing the data for subsequent correla- 
tion calculations, the Fock matrix is diagonalized using the routine DSPEV 
from the LAPACK/BLAS library, upon completion of the self-consistency iter- 
ations. Because the entire spectrum is obtained only once the self-consistency 
has been achieved, it does not strain the computational resources too much. 
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3.21 Subroutine SCF_ UHF 



This subroutine is exactly the same in its logic and structure as the previously 
described SCF_RHF, except that the task of this routine is to solve the UHF 
equation for the system under consideration. Different Fock matrices for the 
up- and the down-spin are constructed and diagonalized in each iteration, un- 
til the self-consistency is achieved. Similar to the case of routine SCF_RHF, 
during the iterations only the occupied eigenvalues and eigenvectors are com- 
puted using the routine DSPEVX. The iterations are stopped once the total 
UHF energy of the system converges to within a user defined threshold. 



3.22 SubrouUne DENSITY 

This subroutine constructs the density matrix for the closed-shell systems, 
assuming orbitals to be doubly occupied. 



3.23 SubrouUne DENSITY_ UHF 

This subroutine constructs density matrices needed for the UHF calculations. 
It generates different density matrices for the orbitals with up (a) and down 
(/?) spins, and in the end adds them to compute the total density matrix. 



3.24 SubrouUne FOCKMAT 

This subroutine constructs the Fock matrix for the closed-shell system and 
calculates the total SCF energy from the density matrix and the one- and 
two-electron integrals. 



3. 25 SubrouUne FO CRM A T_ UHF 

This subroutine is analogous to the routine FOCKMAT, the only difference 
is that it constructs the Fock matrix for the UHF calculations by computing 
separate Fock matrices for up- and down-spins. Thereafter, it calculates the 
total SCF energy from the density matrices and the one- and two-electron 
integrals. 
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3.26 Subroutine ORBDEN 



This routine computes the charge density of a given set of orbitals on user 
specified sites. This helps in understanding as to how the charge is distributed 
in an extended system and helps us to visualize that whether a given orbital 
is localized or delocalized. This analysis can be performed both on RHF and 
UHF orbitals. 



3.27 Subroutine WRITORB 

The purpose of this subroutine is to write the converged orbitals to the orbital 
file in the ASCII format. In case of RHF calculations, orbitals will be written 
in 'ORB001.DAT' file, while in case of UHF calculations, orbitals with up- 
and down-spins will be written in the files 'ORBALPHA.DAT' and 'ORB- 
BETA.DAT', respectively. 



3.28 Subroutine SCI 

This subroutine is the master routine for performing the single configuration 
interactions (SCI) calculations, using the Hartree Fock occupied and virtual 
MOs. SCI is a simple approach which accounts for the electron-correlation ef- 
fects by including configurations which are singly-excited with respect to the 
HF reference state, in the CI expansion[2]. Because of the Brillouin theorem, 
which forbids the mixing of the HF state with the singly-excited configura- 
tions, the ground state remains unchanged in a SCI calculation [2]. However, 
excited configurations do mix with each other leading to a reasonably good 
description of the excited states in this approach[2]. In the present version of 
code, SCI method is implemented only for the case of closed-shell systems, for 
which RHF calculations are performed prior to the SCI ones. To obtain the 
CI eigenvalues, SCI matrix is diagonalized using the routine DSPEV from the 
LAPACK/BLAS library. 



3. 29 Subroutine HAM_ SCI 

This routine constructs the single-CI Hamiltonian matrix for the singlet sub- 
space. A closed-shell reference state is assumed. The Hamiltonian matrix ele- 
ments between different sets of configurations are computed using the standard 
Slater-Condon rules[2]. 
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3.30 Subroutine OPTICS 



This is the master routine meant for computing the optical absorption spec- 
trum, which can be obtained at the single particle level (TB or HF) as well 
as at the SCI level. The range of frequencies over which the spectrum is to be 
computed, along with the line width, are read from the input file. The calcu- 
lated spectrum is written in an output file called 'spec001.dat', which can be 
plotted using xmgrace[34] or gnuplot[35]. At present this routine is restricted 
to the closed-shell systems. 



3.31 Subroutine SPCTRM_ A 



This routine computes the absorption spectrum from the ground state under 
electric-dipole approximation, assuming a Lorentzian line shape and a constant 
line width for all the levels. The computed spectrum is written in an ASCII 
file 'spec001.dat', which can be readily used for plotting using programs such 
as xmgrace[34] or gnuplot[35]. 



3. 32 Subroutine SPCTRM_ lEX 



This is an important subroutine which calculates the linear optical absorption 
of the system from a singly-excited state (as compared to the ground state) for 
one electron TB or HF calculation. The formalism as well as the approach of 
computation is same as that in routine SPCTRM_A, and the output of this 
routine is also written in the ASCII file 'spec001.dat'. Excited state absorption 
spectra are useful in interpreting photo-induced absorption experimental data. 



3. 33 Subroutine DIPMA T 



If optical absorption calculations are to be performed within the SCI approach, 

one needs electric dipole matrix elements between various singly excited con- 
figurations. The task of this routine is to compute these matrix elements so 
that they can be used for calculating the optical absorption spectrum in the 
routine SPCTRM A. 
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3.34 Subroutine NLO 

This subroutine generates the input file to be used for the nonlinear opti- 
cal (NLO) susceptibility calculations at the single particle level. In order to 
reduce the size of dipole matrix elements, one can delete the outermost occu- 
pied and virtual orbitals choosing the card 'ORBDEL'. The program deletes 
the orbitals in the given range using a subroutine named 'delo_nlo', which 
will be described next. Output data of NLO calculations is written in file 
'NLO001.DAT', which includes total number of orbitals, number of occupied 
orbitals, energies of the occupied, as well as virtual states, and the dipole 
moments. NLO option is also restricted to the RHF calculations only. This 
data can be used to compute NLO susceptibilities such as two-photon absorp- 
tion (TPA), third harmonic generation (THG), etc., using the sum-over-states 
(SOS) formalism of Orr and Ward[36]. 

3. 35 Subroutine DEL 0_ NL 

The task of this subroutine is to delete occupied and virtual orbitals from the 
orbital array for the NLO calculations. This routine is called by subroutine 
NLO, if needed. 

3.36 Subroutine DIPCAL 

This routine computes the dipole moment matrix between the single particle 
orbitals under the Hiickel approximation, according to which, among the AOs 
only diagonal matrix elements are nonzero, and their values are nothing but 
their Cartesian locations. This routine simply transforms the dipole matrix 
elements from the AO representation to the MO one. This routine is called 
from the routine NLO and its aim is to supply dipole integrals needed for 
subsequent SOS calculations of NLO susceptibilities mentioned above. 

3.37 Subroutine CI DRV 

This is the driver routine which controls the preparations of various data files 
needed for the post-HF correlated calculations, such as FCI, QCI, MRSDCI, 
etc. The tasks of this routine include: 

(1) Reading the orbitals 

(2) Freezing and/or deleting orbitals if specified by the user in the input file 
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(3) Creating the one-electron effective core potential (ECP) corresponding 
to the frozen set of orbitals, if needed 

(4) Transforming the one- and two-electron integrals from the site (AO) rep- 
resentation to the basis spanned by these orbitals (MOs). They are writ- 
ten to different files so that they can be read and used by a separate 
post-HF correlation program. 



3.38 Subroutine ECP 



This program calculates the effective one-electron frozen core potential matrix 
elements, and stores them in the array. Eventually, these matrix elements are 
added to the other one-electron part of the Hamiltonian before being written 
out in the CI_DRV output files. The energy contribution of the frozen orbitals 
is stored in the variable ecore. 



3.39 Subroutine TWOIND 



This routine performs the two-index transformation on one-electron matrix 
elements, so as to transform them from the AO to the MO representation 
for the use in subsequent correlation calculations. It is called by the master 
routine CI DRV. 



3.40 Subroutine WRITE_1 

This routine writes all the nonzero one-electron Hamiltonian matrix elements 
expressed in the MO representation in a binary file called 'ONEINT001.DAT'. 
This routine is also called by the master routine CI_DRV, and the output is 
meant to be used in the post-HF correlation calculations. 



3.41 Subroutine FOURIND 

This routine performs the four-index transformation on the two-electron in- 
tegrals so as to transform them from the AO to the MO represention. This 
routine is also called by the master routine CI DRV. 
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3.42 Subroutine WRITE_2 



This routine writes all the nonzero two-electron Hamiltonian matrix elements 
expressed in the MO representation in a binary file called 'TWOINT001.DAT'. 
This routine is also called by the master routine CI_DRV, and the output is 
meant to be used in the post-HF correlation calculations. 

343 Subroutine DIPINT 

This routine transforms the dipole operator from the site (AO) representa- 
tion into the representation of the orbitals (MOs). The routine is called from 
CI_DRV and it provides dipole matrix elements to be used for optical prop- 
erties calculations in conjunction with post-HF correlated calculations. The 
algorithm used in this routine is the one described in the routine DIPCAL. 

344 Subroutine DIP OUT 

The purpose of this routine is to write all the nonzero dipole matrix elements in 
the binary file 'DIPINT001.DAT' to be used in conjunction with the post-HF 
correlated calculations. This routine is also called from CI DRV. 



4 Installation, input files, output files 

We believe that the installation and execution of the program, as well as prepa- 
ration of suitable input files is fairly straightforward. Therefore, we will not 
discuss these topics in detail here. Instead, we refer the reader to the README 
file for details related to the installation and execution of the program. Ad- 
ditionally, the file 'manual.pdf' explains how to prepare a sample input file. 
Several sample input and output files corresponding to various example runs 
are also provided with the package. 



5 Results and Discussions 

In this section, we present and discuss the numerical applications of our code. 
First we present the results on the convergence of total energy with respect 
to the oligomer's length for various conjugated polymers such as i-PA, PPP, 
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PPV, and polyacene. We also provide our results for the UHF calculation on 
graphene nanodisks. Next, we present our results for the tight binding band 
structure calculations for PPP. Further, we discuss and compare the linear 
absorption optical spectrum of PPP computed by different approaches such 
as TB, HF, and SCI. As mentioned in the above sections that the program can 
perform calculations in the presence of finite electric field as well, which makes 
it capable of calculating electro-absorption spectrum. For centro-symmetric 
systems, the electro-absorption spectrum allows one to explore both the even 
and the odd parity excited states in a linear optical absorption process. Results 
and analysis of the electro-absorption spectrum of an oligomer of PPP at the 
SCI level are also presented. 

5.1 Total Energy 

Our code can be used to study both the ground and the excited state properties 
of oligomers of various polymers because they are nothing but finite molecules, 
ranging in size from small to large. However, in this section we demonstrate 
that our code can also be used to obtain the ground state energy/cell, in the 
bulk limit, for one-dimensional periodic systems such as polymers. 

The energy per unit cell of a one-dimensional periodic system can be obtained 
using the formula 

Eceii = A^n = - E^), (13) 

where En+i/En represent the total energies of oligomers containing n + 1/n 
unit cells. Thus, using this formula, for sufficiently large value of n, one can 
obtain the energy/cell of a polymer in the bulk limit, from oligomer based 
calculations. In what follows we show that value of Eceii converges quite rapidly 
with respect to n. 

First we examine the convergence of Eceii obtained using Eq. (13) with respect 
to the number of unit cells n, using the RHF approach. The energy convergence 
threshold for the SCF calculations is set up to eighth decimal place. In Fig. 1 
we plot A.En as a function of n, for various polymers such as t-PA, acene, PPP, 
and PPV. From plots (c) and (d) in Fig. 1, it is evident that energy for PPP 
and PPV converges rapidly and convergence is achieved for short oligomer 
length (6-7 units), while in t-PA and acene (Fig. 1(a) and (b), respectively), 
the convergence is achieved for larger value of n. 

The above mentioned calculations can also be performed using the UHF ap- 
proach. However, it is more meaningful to use this approach for systems with 
open-shell configurations, such as trigonal zigzag graphene-nanodisks. Thus, 
UHF calculations can be used to ascertain whether or not a high-spin state is 
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Figure 1. Convergence in AE^ with respect to the ohgomer length (n) of various 
systems: (a) t-FA, (b) polyacene, (c) PPP, and (d) PPV. Inset of each graph shows 
the zoomed plot, in the energy window of lO"'^ - 10~* eV. 



Figure 2. Schematic diagram of the trigonal zigzag graphene nanodisk with six ben- 
zene rings (22 sites) considered in this work. 

the ground state. For example, we found that the ground state of a trigonal 
zigzag graphene-nanodisk with six benzene rings (c/. Fig. 2) is a triplet one 
(see Table 1). 



5.2 Band Structure 



Our code can be used to perform band structure calculations using the nearest- 
neighbor TB model. In Fig. 3, we present band structure of PPP. It possesses 
six bands out of which two are localized (labeled as 'Z', 7*'), while four are 
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Table 1 

Results of UHF total energy calculations on a trigonal zigzag graphene nanodisk 
with six benzene rings (c/. Fig. 2) for various combinations of majority (no-) and 
minority (n^j) spins. Our results indicate that the triplet state being significantly 
lower in the total energy than the singlet state, is the true ground state. 
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Figure 3. Band structure of PPP, computed using Hiickel model Hamiltonian. Dis- 
persion less flat bands denote the localized levels (l/l*), while the other bands rep- 
resent the delocalized ones. 



delocalized (labeled as ^di/d2\ 'rfi/rfg')- As is obvious from the Fig. 3, the 
band structure has a direct band gap, Eg = 2.16 eV at the F-point. Our 
calculated band structure is fully consistent v^^ith that reported by Heeger and 
co-workers [37]. 
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5.3 Linear optical absorption spectrum 



In order to understand the optical properties of conjugated polymers, it is 
essential to study their low-lying excited states. The linear optical absorp- 
tion spectrum, which is all about one-photon absorption, provides an efficient 
medium to explore the odd-parity low-lying states in centro-symmetric sys- 
tems. In an independent particle model it can be explained in terms of excita- 
tions from ground state, by promoting one electron from one of the occupied 
MO (valence band states) to one of the unoccupied MO (conduction band 
states). While, in a correlated electron approach, such as the P-P-P model, 
this one-electron picture does not hold, and the excited states are linear com- 
bination of several configurations obtained by performing CI calculations. 

Our code can be used to compute linear optical absorption spectrum of tt- 
conjugated materials using the Hiickel model (TB approach), as well as P-P-P 
model Hamiltonians employing both the HF and the SCI approaches. 

In Fig. 4 we present the linear optical absorption spectrum of planar PPP with 
eight repeat units (PPP-8) using TB, HF, and SCI approaches. The spectrum 
was calculated using the carbon-carbon bond length (within a phenyl ring) 
to be 1.4 A, while the inter-ring single bond length was taken to be 1.54 A. 
The standard value of -2.4 eV was used for the intra-benzene hopping, while 
for the inter-ring single bond we used the value -2.23 eV. The HF and the 
SCI calculations were performed using the standard parameters in the P-P-P 
Hamiltonian. All the three graphs depicted in Fig. 4 are qualitatively quite 
similar, with the first two features, and the most intense feature, corresponding 
to the x-polarized transition. The spectrum drawn at the HF level using P-P-P 
model is blue shifted as compared to the Hiickel model spectrum, which is due 
to the well known tendency of the HF approach to over estimate the band gaps. 
However, inclusion of correlation effects at the SCI level brings the spectrum 
back to lower energy range. The location of the lowest peak corresponds to the 
optical gap, which is found to be around 3.5 eV, in excellent agreement with 
the reported first-principles calculations [38], and also in good agreement with 
the experimental value (~3 eV) for the infinite polymer [39]. A higher level CI 
calculations like FCI, QCI, and MRSDCI can also be performed using, e.g., 
the MELDF code [40] beyond HF. The input for these higher CI calculations 
can be prepared using CIPREP option in the input file. 

On analyzing the wave functions of the excited states contributing to the 
spectra, we find that first peak in TB and HF spectra is due to transition 
from the highest occupied molecular orbital (HOMO) to the lowest unoccupied 
molecular orbital (LUMO), i.e., H ^ L. In the spectrum obtained using SCI 
approach, for peak I (see Fig. 4), the major contribution comes from H ^ L 
configuration, while the contributions due toH— l^L + 1, H — 2^L + 2, 
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Figure 4. Linear optical absorption spectrum of planar PPP-8, computed using the 
Hiickel as well as P-P-P models, both at the HF and the SCI levels. Various peaks 
and their polarizations are labeled. The spectra are plotted using the line width of 



H — 3 —>^ L + 3, and H — 12 L + 12 are also significant. 

Our code is also capable of examining the role of Coulomb parameters used in 
the P-P-P Hamiltonian by not only considering standard and screened set of 
parameters defined earlier (c/. Sec. 2), but also any user defined set of Coulomb 
parameters, which may be material specific. Fig. 5 depicts the linear absorp- 
tion spectra computed using standard parameters (dashed line) and screened 
parameters (solid line), for planar geometry of PPP-8 at SCI level. Both the 
spectra are qualitatively similar, while quantitatively, spectrum computed us- 
ing screened parameters is red shifted as compared to the standard parameter 
spectrum. We note that the value 3.27 eV of the optical gap obtained using 
the screened parameters is in better agreement with the experimental value 
(~ 3eV) for the infinite system[39]. 

5.4 Electro-absorption spectrum 

As mentioned earlier our code can also solve the HF equations in the pres- 
ence of a finite electric field. This allows us, in principle, to use it for com- 
puting quantities such as the static dielectric polarizabilities and the electro- 
absorption (EA) spectra of various conjugated systems. For example, the first- 
order polarizability is determined as the second derivative of the total energy 



0.1 eV. 
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E(eV) 

Figure 5. Linear optical absorption spectrum of PPP-8 computed using the standard 

(dashed Hne) as well as the screened (solid line) set of Coulomb parameters in the 
P-P-P model, using the SCI approach. The spectra are plotted using the line width 
of 0.1 eV. 

per unit cell (Etotai) with respect to the strength of the external electric field E. 
The second derivative of £'totaZ can be calculated employing finite difference 
formulas, using Efotai with, and without, the electric field. However, in this 
work we are concentrating on another possible application of the finite-field 
approach, namely the electro-absorption spectrum. 

The electro-absorption spectrum is an important experimental probe for var- 
ious materials which, through linear optical absorption allows, one to probe 
both even, as well as odd parity states for centro-symmetric systems. The note- 
worthy point is that without electro-absorption normally one will have to use 
both linear absorption as well as nonlinear absorption to investigate states of 
both parities. Using our program one can easily make theoretical prediction of 
this quantity using its definition: Aa{uj, E) = a{u!, E) — 0(0;, 0), where a{uj, E) 
is the optical absorption spectrum in the presence of the external electric field 
E. In order to compute it, the HF equations are solved in the presence of the 
field E (c/. Eq. 6) leading to a modified set of orbitals and their eigenvalues, 
as compared to the zero-field (E = 0) case. This orbital set can subsequently 
be used to compute: (a) optical absorption at the HF level, or (b) optical 
absorption at any level of correlation. Then by using the above mentioned 
definition of Aa(cj,E) one can compute the corresponding EA spectrum. In 
this work we present the results of the EA spectrum of PPP-8, computed at 
the SCI level using our program. The EA spectra are usually interpreted in 
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Figure 6. Left panel: (a) The linear absorption spectrum planar molecule of PPP-8, 
computed in the presence of finite electric field, using standard parameters in the 
P-P-P Hamiltonian with the SCI approach, (b) Electro-absorption spectrum (solid 
line) and the first derivative of linear absorption with respect to the photon energy 
E (dashed line) for the same system. Right panel: The two photon absorption (TPA) 
spectrum of PPP-8, computed at SCI level, using the MELDF program[40]. In all the 
figures above, E represents the photon energy. The mAg state in the TPA spectrum 
corresponds to the encircled feature in the EA spectrum, as is also clear by the peak 
positions. 

terms of Stark shift analysis and it roughly follows the first derivative of linear 
absorption, calculated in the presence of the electric field, with respect to the 
photon energy. 

Fig. 6 (a) shows the linear absorption spectrum computed in the presence of 
finite electric field, using the standard parameters in the P-P-P Hamiltonian 
with the SCI approach. As expected, the most intense features in the EA 
spectrum appears when the field is applied along the conjugation direction 
which happens to be the x-axis in the present case. In these calculations 
the electric field strength was taken to be 10~^ V/A, which is of the same 
order of magnitude used in various experiments[41]. Fig. 6 (b) depicts the EA 
spectrum together with the first derivative of linear absorption with respect 
to the photon energy E. We note that the EA spectrum resembles the first 
derivative of linear absorption spectrum, roughly. The first absorption peak 
in linear absorption spectrum occurs at 3.56 eV, following the other features 
at 4.21 eV, 5.8 eV, and 6.27 eV [see Fig. 6 (a)]. In EA spectrum, the first 
peak appears at 3.51 eV, which is similar to da/dE spectrum [dashed curve in 
Fig. 6 (b)]. We therefore assign this feature to a Stark shift of the lowest odd 
parity (l-B^) exciton. There is also a feature at 6.09 eV [circled in Fig. 6 (b)], 
which cannot be found by a derivative analysis of the absorption spectrum, 
and must therefore be due to an even parity (mAg) exciton. To verify this, 
we also calculated two photon absorption (TPA) spectrum using the MELDF 
code[40] at the SCI level. The right panel of Fig. 6 shows our calculated TPA 
spectrum for PPP-8. The spectrum features an intense peak located at 3.04 
eV, which exhibits strong dipole coupling with the lowest one-photon IBu 
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state, and therefore is the mAg state[42], located at 6.08 eV. Thus, our result 
for mAg state obtained from the EA spectrum is in excellent agreement with 
that from the TPA. 



6 Conclusions and Future Directions 

In this paper we have described our Fortran 90 program which solves the 
HF equations for both the closed- and open-shell molecular systems using the 
semi-empirical Hiickel and PPP models. To demonstrate the features of our 
code, we presented results of numerous test calculations on various molecular 
systems. The reason behind developing the present program is twofold: (a) to 
develop a code in a modern language such as Fortran 90 which can carry out 
dynamic array allocation, and thus, free the user from specifying and changing 
array sizes, and (b) to provide an open software which will be widely available 
to users which they can use and modify as per their needs. One possible 
generalization of this code will be to include the capabilities for calculations on 
infinite systems in one- and higher dimensions to calculate the band structure 
and related properties using the P-P-P model. Work along those directions 
is continuing in our group, and results will be published as and when they 
become available. 
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